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Abstract 

A lattice version of the Fokker-Planck equation (FPE), accounting for dissipative interactions, 
not resolved on the molecular scale, is introduced. The lattice FPE is applied to the study of 
electrorheological transport of a one-dimensional charged fluid, and found to yield quantitative 
agreement with a recent analytical solution. Future extensions, including inelastic ion-ion collisions, 
are also outlined. 
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I. INTRODUCTION 



Over the last decade discrete lattice versions of kinetic equations, most notably the 
Lattice-Boltzmann (LB) method, have undergone burgeoning progress for the simulation 

nnnn 

of large scale hydrodynamic flows (H, [2j, y, U] and of the dynamics of colloidal suspensions 



0,3. 



The LB method was developed in response to major flaws of its predecessor, the 



lattice gas cellular automaton |2Ll£(|, following the canonical route of classical Statistical Me- 
chanics. A few years later the formal connection of the LB method with continuum kinetic 
theory was also elucidated 3. 

One of the major appeals of the LB method is its flexibility, which allows to accomodate a 
host of complex physical effects, including boundary conditions at interfaces, intermolecular 
forces and even chemical reactions, through efficient and elegant discretizations of the force 
term in the kinetic Vlasov-Boltzmann equation. One of the limitations of the LB method 
is that, since it generates the time evolution of the one particle distribution function, fluc- 
tuations are not accounted for. "Brownian noise" becomes increasingly important as one 
explores flows on ever smaller scales, as in colloidal systems, or in narrow pores (microflu- 
idics). The problem has been addressed on the colloidal scales, by incorporating a random 
(Brownian) force component in the momentum flux stress tensor , or at the level of 

the discrete LB equation itself 

On the other hand, all LB implementations to date are based on the relaxation form of 
the co Uisi „„ operas, either in ^ Q 0r tensoria! ta Q. In this paper we show that 
the LB methodology can be easily extended to deal with small scale processes involving a 
Brownian component, by introducing a lattice version of the Fokker-Planck (FP) collision 
operator Q| In full analogy with the LB equation, the present lattice Fokker-Planck 
equation builds upon an optimized form of importance sampling of velocity space which, at 
variance with most numerical grid methods [lj, Q Q|- permits to solve the Fokker-Planck 



equation near local-equilibrium in full single-particle phase space. The FP operator accounts 
for frictional dissipation due, e.g., to solute-solvent interactions or inelastic collisions with 
obstacles and confining surfaces, which need not be resolved on a molecular scale. We 
have, in particular, in mind electrorheological transport of ions confined in swollen clays, 
between membranes, or through water-filled nanopores (see e.g. jjjj]), like ion channels 
through membranes The FP equation has recently been applied to a detailed analysis 
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of stationary ion currents through one-dimensional pores of finite or infinite length [121, and 
analytic results were obtained in the case of independent (non-interacting) ions. In particular 
the ion current was shown to saturate with increasing applied field, and this result will serve 
as a benchmark for the numerical calculations presented later in this paper. 

The paper is organized as follows. In section 2 we discuss the respective roles of the FP 
and Bhatnagar-Gross-Krook (BGK) collision operators in modelling transport and flows of 
solutes in implicit solvent and confined geometries. A one-dimensional kinetic model for ion 
transport through narrow pores is presented in Section 3. The lattice Fokker-Planck (LFP) 
equation for this problem is derived in Section 4. Numerical results are presented in Section 
5 and compared to the predictions of ref. jl^, while concluding remarks are made in Section 
6. Stability analysis of the LFPE is outlined in the appendix. 



II. FOKKER-PLANCK VERSUS BGK COLLISION OPERATORS 

Kinetic equations for the time evolution of the distribution function /(r, v; t) of a given 
species of particles conventionally involve a free-flow term (left-hand side) and a collision 
term (right-hand side), i.e. 

(d t + v a d Xa + a a d Va )f(r, v; t) = C[f(r, v; t)} (1) 

where x a and v a (1 < a < d) are the cartesian components of the d-dimensional position 
and velocity vectors r and v, d Xa and d Va are the components of the corresponding gradient 
operators, while a a = F a /m are the components of the acceleration due to an externally 
applied, or a self-consistent force field F (m is the particle mass). The Einstein convention of 
summation over repeated indices a is assumed. On the right hand side C denotes a collision 
operator (yet to be specified), acting on the distribution function /. 

If one is interested in flows involving two (or more) species, e.g. a solvent and a solute, 
one may adopt one of two strategies: 

a) one may treat the two species (say a and b) on the same footing, by introducing two 
distribution functions fi(r, v; t) (i = a or b) , associated with the two species, the time 
evolution of which is governed by two coupled kinetic equations. In the perspective of a 
discrete lattice formulation of the LB type, the simplest collision operators Cij(i,j = a or b) 
acting on the distribution functions are of the BGK form, involving several relaxation times 
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m 

b) An alternative route, which is particularly appropriate when the solutes are colloidal 
particles, and which will be adopted here, is to assume a separation of time scales, and use 
an implicit solvent description of the Langevin form, involving frictional and random forces. 

nn 

The corresponding collision operator is the familiar Fokker-Planck operator |13L Il4| which 
may also account for inelastic collisions of the solute particles with confining surfaces. In 
such an effective one-component description, the dissipative solute/solvent and solute/wall 
couplings are described by a FP collision operator, while conservative solute/solute collisions 
are described by a BGK operator. Explicitly, in the presence of the electric field, 

(d t + v a d Xa + a a d Va )f(r, v; t) = C FP [f(r, v; t)} + C BGK [f(r, v; t)} (2) 

where 

C FP [f}=d Va (R a + Dd Va )f (3) 

C BGK [f] = -u(f - f B % K ) (4) 

In eq. © R a are the components of the drag force, which will be taken of the familiar 
form R = TV, where t is a constant friction coefficient, D characterizes diffusion in velocity 



space, and is related to t by D = jv^, where vt = ykBT/m is the thermal velocity and 
uj = 1/r is the solute-solute collision frequency. 

The BGK collision operator entails a relaxation of the distribution function to the usual 
local Maxwellian equilibrium, 

J BGK\ i i ) (2^2)3/2 e ^ 

where n(r; t) is the local density, u = u(r; t) is the local average (or flow) velocity, while the 
thermal velocity vt may depend on position and time, via the local temperature T(r;t). 

The equilibrium resulting from the combined action of Fokker-Planck collisions and the 
electric field is given by a global shifted Maxwellian: 

Wlr,v,tj- (27ru 2 )3/2 e [b) 

where u E = ^ is the drift speed associated with the electric field. By definition, the com- 
posite local equilibrium resulting from the competition of the two operators in the presence 
of an electric field satisfies the condition 

c fp \n + c bgk m - a a o Va r = o (7) 



This may be rewritten as 

c FP ' E [n = ^(r-r B q GK ) (s) 

where C FP ' E = d Va (v a + v^d Va ) — (a a /j)d Va is the Fokker-Planck operator including the 
electric field. 

From the expression (JBJ), it is clear that in the strongly dissipative regime cu/7 — > 0, 
f eq ~* fi?PEi while the leading order correction is proportional to the difference fpp e — fBGK 5 
i.e. to the deviation of the local flow field u from the drift speed Ug. This is consistent with 
the fact that when u = u^, a condition which is reached at steady-state, the FP and BGK 
equilibria coincide. 

The interplay between the FP and BGK collision operators spawns a rich variety of 
physical effects, which will be the object of future investigations. In the sequel, however, we 
shall confine our attention to the methodological aspects related to the lattice formulation 
of the Fokker-Planck equation. For illustration purposes, and the sake of simplicity, we 
henceforth restrict the discussion to the case d = 1 (one spatial dimension) of the FP 
equation and to its validation by comparison with recent exact results. 



III. A ONE-DIMENSIONAL KINETIC MODEL FOR TRANSPORT THROUGH 
PORES 

In this Section the implicit solvent kinetic model introduced in the previous Section is 
applied to the important problem of single-file ion transport through a water-filled pore 
connecting two reservoirs, under the action of an applied electric field or ion-concentration 
gradient. The one-dimensional version of the kinetic model ©-Q provides a crude repre- 
sentation of ion permation of ion channels through membranes separating intra and extra- 

nn 

cellular compartments |lja, Ion permeation of such channels has been examined by 

numerous Molecular Dynamics (MD) or Brownian Dynami cs ( BP) simulations of realistic 
or semi-realistic quasi-cylindrical models (for a review, see (2l|) or by numerical solutions 
of the Poisson-Nernst-Planck equations |22fl, but the present study is inspired by the recent 
kinetic modelling of ref. 

The action of the confining, quasi-cylindrical pore is crudely represented by restricting 
ion motion to one dimension and by a contribution to the frictional force —^yv. 
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The one-dimensional version of the kinetic equation ©-(III) may be cast in the form 

(d t + vd x )f = C FP [f] + C BGK [f] - adj = dMvf + v 2 T d v f)} - u(f - f« GK ) - adj (9) 

where / = f(x,v;t) and a = qE/m, E being the applied electric field and q the charge of 
the ions; in practice one is mostly interested in mono or divalent cations (q = +e or +2e). 
The last two terms on the r.h.s. of eq. (JHJ) may be regrouped into a single FP-like term 

C FP [f] - adj = d vl [(v -u E )f + v 2 T d v f] (10) 

where ue = qE/m'j = a/7 is the ion drift velocity in response to the applied field. 
The local equilibrium solution of the kinetic equation (J^J is hence 

JBGK\ X i V i ^) ~ ( 27lv 2 T y/2 e ' ' 

The zeroth, first and second moments of the distribution are the local density n, current J 
and pressure (or momentum flux) P per unit mass 



n(x;t) = / f(x,v;t)dv (12) 

/oo 
vf(x,v;t)dv = n(x;t)u(x;t) (13) 
-00 

POO 

P(x;t) = / vvf(x,v;t)dv (14) 



Note that in one dimension the momentum flux is proportional to the (kinetic) energy, but 
this is of course no longer true in higher dimension. 

By multiplying both sides of the kinetic equation Q successively by 1, v and v 2 , and 
integrating over all v, one easily arrives at the following macroscopic equations 

d t n{x;t) + d x J{x;t) =0 (15) 

d t J(x; t) + d x P(x; t) = -«fJ(x\ t) + n(x; t)a (16) 

d t P(x; t) + d x Q(x; t) = -2 7 n(x; t) + 2a J{x; t) (17) 

2 

Eq. lfTBIl is the continuity equation expressing the conservation of mass; equation (fTfil) ex- 
presses momentum balance with proper account of friction and acceleration due to the 
electric field while (fT?]) is the energy balance equation, taking into account the heat flux 

/oo 
v 2 vf{x 1 v]t)dv (18) 
-00 



as well as frictional dissipation. The standard kinetic equation with the BGK collision 
operator yields the same three macroscopic equations without the dissipative contributions 
stemming from the FP collision operator. In the case of a steady, homogeneous flow, eq. (fT6l) 
leads back to Ohm's law, qJ = aE, with a conductivity a = nq 2 /7m. 

Stationary solutions of eq. (jUJ) have been obtained in the case of independent ions (00 — 0), 
for finite length as well as infinitely long channels in ref. Solutions in the non-stationary 
case, and including the BGK collision term in addition to the FP term, can only be obtained 
numerically. In the next section we derive the lattice version of the kinetic equation (JHJ). 



IV. THE LATTICE-FOKKER-PLANCK EQUATION 



The lattice version of the kinetic equation Q may be systematically derived along the 

Q finnn 
. Since the latter is well documented[l|,|2j,|3j,|4| we restrict 

most of the following discussion to the FP collision operator (fT0|) . 

The distribution function is expanded onto a Hermite basis 



K 



f(x,v;t) = J2 F k(%;t)h k (v)w(v) 

k=0 



(19) 



where w(v) = {2Txv 2 T )- 1 ' 2 e- v2 / 2v T is the one-dimensional Hermite weight function, while 



hk(v) is the Hermite polynomial of order k. By substituting eq. (fT9|) into the kinetic equation 
eq.©, and projecting upon the Hermite basis, one arrives at 

d^x; t) + d^x; t) = Q(x; t) (20) 

where 

/oo 
f{x,V]t)hi(v)dv (21) 
-00 

/oo 
vf(x,v,t)hi(v)dv (22) 
-00 

/oo 
C FP [f(x,v;t)Mv)dv (23) 
-00 

C FP being the linear FP operator (jlOj) . 

Equations (J2Sj) are the usual moment relations associated with the FP equations. The 
next step is to evaluate the kinetic moments Fi, Gi and Ci by Gauss-Hermite quadrature. 
Noting that f(x,v;t)/w(v) is a polynomial in v, the quadrature reads 



Fi(x;t) = ^—j-^ w ihi( v i) 

i=0 w \ v i) 



(24) 



where the Vi and Wi are the nodes and weights of the quadrature. It is crucial to observe 
that Eq. is exact for polynomials of degree up to (2G + 1). This means that the 

present Hermite-Gauss projection is de-facto equivalent to an optimized form of importance 
sampling of velocity space. It is this discretization which permits to compute the evolution 



of the low-order macroscopic moments 
discretization of velocity space (see ref. 



density-current) more efficiently than any standard 
15| and references therein). In fact, any grid method 
would necessarily use a mesh-spacing in velocity space which is a fraction, say at most 1/10, 
of the thermal speed, vt- Covering a few units of vt would then take at least 30-50 grid 
points in velocity space. In contrast, LFPE only needs at most five. The advantage to LFPE 
would be even more substantial in higher dimensions. 

Expressions similar to Eq.""""[|) hold for Gi and C\. Substituting these into eq. ("2*njL and 
identifying the factors of hi(vi) on both sides of the resulting sum, one obtains the following 
set of equations 

d t f i (x;t)+v i d x f i (x;t) = a(x;t) < % < G - 1 (25) 

where the following identifications have been made 

^■•t^Htir (26) 

C , (I;( ),Q%* (27) 

W[Vi) ' 

The discrete collision operator is entirely specified by the coefficients q defined by eq. (""T"") . 
These can be unambiguously computed from the spectral decomposition of the continuous 
operator C FP [f(x, v; t)] similar to eq. ("H""), i.e. 

C FP { Vi ) = C FP [f(x,v t -t)] = J2Ck(x-t)h k ( Vl )w(v t ) (28) 

k=0 

Knowledge of the spectral coefficients Ck(x,t) allows the discrete coefficients Ci(x;t) in 
eq. ((2*7)1 to be calculated, thereby providing an operational definition of the discrete FP 
operator. 

The last step is to perform time integration along the characteristics dxi = Vidt according 
to the standard practice of LB algorithms. This yields the desired lattice-Fokker-Planck 
(LFP) equation 

fi(x + Vi At, t + At)- fi(x, t) = a(x; t)At (29) 
8 



where At is the time-step chosen for the numerical solution. 

The spectral coefficients of the FP operator are easily calculated to be 



C = (30) 

C x = -7 J + na (31) 

C 2 = -2 7 (P - nvl) + 2aJ (32) 

C 3 = -&y{Q - J) + 3a(P - nu£) (33) 

C A = -4 7 (P - 5^P + 2ra4) + 4a(Q - 2Jv%) (34) 

where P = $™ 00 dvv A f(x 1 v]t) and the first four Hermite polynomials are h = 1, hi = v, 



h 2 = v 2 — h 3 = v 3 — 3ffy, hi = v 4 — 4u 2 t>y + t^- High order coefficients fc > 2) 
differ from those derived from the BGK operator (if uj = 7). All coefficients, as well as n, 
it, P and Q, defined by eqs (|T2"|) - (|TH) and (jTHjl . and P, are local quantities, depending on x 
and £. 

Returning to the LFP equation (|2H . we now consider two discretized models, namely 
Q 3 and Q 5 , involving three and five velocities and where v\ = 1/3 and v% — 1 respectively. 
Within Q 3 the velocities t>j and associated weights Wi are 

f = wo = 2/3 
U! = +l; wi = 1/6 (35) 
v 2 = -1; u> 2 = I/ 6 

The three velocities t>j (i = 0, 1,2) may be considered as the components of a three-vector 
(0, 1, —1). It proves convenient to construct a basis of three such vectors which satisfy 
the orthonormality conditions 

A« ■ A« = f>f = (36) 

i=0 

These vectors are easily calculated to be 
A<°> = (1,1,1) 

A« = V3(v , Vl , v 2 ) = \/3(0, 1, -1) (37) 

a (2) = ^{vi-viy x -viy 2 -vi) = l={-w) 
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within this basis, the collision operator in eq. (|29j) may be cast in the form 

d(x, t) = (Co4 (0) + + C 2 Af V< (38) 



The various x and t-dependent quantities follow from the definitions (|12| - (fT4j) by replacing 
the integrals over velocities by sums over the 3 discrete velocities t>o, V\ and v 2 , using the 
proper weights u>j. These prescriptions entirely specify the Lattice Fokker-Planck (LFP) 
equation for the model Q3. The same procedure may be used to specify the LFP equation 
for the Q$ model, where the 5 velocities are (—2, —1, 0, 1, 2), and the corresponding weights 
Wi are (1/12, 1/6, 1/2, 1/6, 1/12), with a resulting thermal speed vt — 1. 

The Lattice BGK equation has been widely described in the literature j^j], and therefore 
only a very sketchy description is provided below. The lattice BGK equation has the same 
form as eq.JSHJ), where the lattice BGK collision operator is given by 

cf GK = -u(fi ~ fsGK,) (39) 

The local equilibrium takes the form of a local Maxwellian expanded to second order in the 
Mach number M = u/vt, 

f bgk,% = w<n\l + ^ + J (40) 

The lattice equilibria, including the weights Wi, are designed on the requirement of fulfilling 
mass and momentum conservation, that is 

i=0 i=0 

and 

G-1 G-1 

Z! f'BGK^i = X M = UU ( 42 ) 
i=0 i=0 

The parameter to is an inverse relaxation scale which controls the shear viscosity of the 
lattice fluid. 



V. NUMERICAL RESULTS 

As a first application and test of the LFP formalism, we consider the single-file trans- 
port of ions through a pore of length L connecting two reservoirs containing ions at given 
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concentrations. As in ref. |3| we assume that at both ends of the pore, the ion distri- 
bution function is a Maxwell-Boltzmann equilibrium distribution at temperature T which 
determines the thermal velocity Vt- For instance, in the Q3 case, boundary conditions are 
imposed as follows 

/,(* = <>,*) = J= 

f 2 (x = L + l,t) = -£= (43) 

V Z7T 

where subscripts 1 and 2 stand for rightward and leftward propagation and n/, n r indicate 
the left and right reservoir densities respectively. The reservoirs are located at x = and 
x = L + 1 respectively, while the physical channel runs from 1 < x < L. Note that the 
expressions (jlHJ) correspond to fixing the incoming fluxes from the reservoirs, hence they 
do not imply that n\ and n r coincide with the fluid density at inlet and outlet sections. In 
fact, since the ion distribution function in the reservoirs is a Maxwellian at zero macroscopic 
speed, continuity of the fluxes (current density), implies a discontinuity of both density and 
velocity profiles at both ends of the channel. 
We define the dimensionless current 

J* = jM- (44) 
V v T n 

the dimensionless acceleration, or electric field 

maL qEL E . . 

a = k^T = -kjr = e; (45 » 



and the dimensionless collision rate 



lL 7 £ 7 



= — = — = # (46) 
Vt 7t Et 

where n is the reservoir ion density, Et = ksT/qL is a "thermal" electric field, such that the 
work it produces to move a charge q over the channel length L equals the thermal energy 
ksl ', while E y = ^ksT / '(qvr) = m^VT/q is the electric field producing on a charge q a force 
which balances the frictional force —m^yvT- In other words, in a linear regime, the drift 
velocities associated with Et and E 1 are ut = ^t/t* and w 7 = vt respectively. Clearly, 
in infinitely long channels, or zero-temperature fluids, Ut — > 0, and the ionic current is 
controlled by pure dissipation, ue — In finite-size, finite-temperature situations, and 



11 



constant-flux boundary conditions, however, deviations from this simple Ohmic regime must 
be expected, as we shall show in the sequel. 

We have solved the LFP equation numerically to determine the stationary distribution 
function f(x,v) and derive the corresponding macroscopic moments, primarily the current 
J(x). Since we focus on the Fokker-Planck operator, for the time being, we exclude ion-ion 
collisions by setting u = 0. 

In Figure 1, the stationary density profiles n(x) calculated with the Q3 and Q$ models 
for two values of the electric field, a* = 1 (moderate) and a* = 3 (strong), are compared to 
the analytical result of the continuous velocity model [3], which reads as follows: 



n eX act(x) = Ae a * x/L + B 

where 



A 



2 sinh(a*/2) 



a 1 — 

n r — n>i H y2ixB 

a* 



nie a * I 2 — n r e a * I 2 

B 



2sinh(a*/2)e- a *W 2 + x /f^ cosh(a*/2) + sinh(a*/2) ZTJ^h* dv ®( v 



and $ = J l/27Tfrexp(— v 2 j2v\). The total current flowing through the channel is J exac t = 
aB/j. It should be appreciated that the strongly non-linear dependence of the current (see 
the expression of B) reflects the highly non-trivial competition between the effect of the 
boundary conditions and the electric field. 

The ratio of the numerical to analytical density profiles n(x) / n exact (x) , is seen to deviate 
from the unit value of about 10% for a* = 3 and less than 5% for a* = 1. These departures 
can be attributed to the constant-flux boundary conditions which, as discussed previously, 
introduce a discontinuity at the open ends of the channel. We also notice that the Q 5 
solution appears appreciably more accurate than the Q3 one. At lower values of the electric 
field, both Q 3 and Q 5 yield excellent agreement with the analytical solution. 

The current-voltage relation is illustrated in Fig. 2, where the reduced current J* is 
plotted as a function of the reduced electric field, a* = E/Et, for three values of the 
reduced collision rate 7* = E 1 /E T = 1,5, 10. Results obtained with the Q3 and Q5 models 
are compared to the analytical predictions of ref. [3], where the problem was solved for a 
continuous velocity spectrum (Qoo)- First, we notice that the Ohmic regime, E/E T <ti 1, 
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is well reproduced in all three cases. As the strength of the electric field is increased, the 
current shows the saturation effect predicted by the analytical results 

This saturation is a direct consequence of the imposed boundary conditions. Since the 
reservoirs cannot feed more than uvt ions per unit area and time, the current cannot exceed 
its ballistic value VTnj^hx in the limit E — > oo. The rate of convergence towards the 
asymptotic limit depends sensitively on the reduced collision rate 7*. Fig. 2 shows that the 
results obtained with the Q$ model are systematically closer to the exact results than those 
obtained with the Q3 model. 

We have also checked that with simple boundary conditions, e.g. periodic, Ohm's law 
J* = E*/j*, is reproduced to machine accuracy and independently on the grid size, beyond 
about 100 grid points. However, the present case is much more challenging due to the highly 
non-trivial boundary conditions. In fact, a very slow convergence of the results as a function 
of grid resolution is observed (see Figure 3). This is probably caused by the representation of 
the reservoirs by a single lattice point without making allowance for any finite size transition 
region between the current-carrying distribution function in the channel and the no-current 
equilibrium current in the reservoirs. The consequence of a less abrupt treatment of the 
boundaries will be explored in the future. 

VI. EFFECTS OF ION-ION COLLISIONS 

The BGK collision operator (|39|) does not make any contribution to the continuity and 
momentum equations because, by construction, ion-ion collisions have been designed to 
conserve mass and momentum (elastic collisions). As a result, the BGK operator does not 
contribute to the momentum equation, and consequently, it cannot produce any effect on 
the electric conductivity of the system. In fact, ion-ion collisions are in control of the shear 
viscosity of the charged fluid, clearly an irrelevant notion in one-dimensional systems. 

However, in electro-rheological flows, ions move in the presence of a surrounding solvent, 
typically water. It is therefore reasonable to assume that ion-ion "effective" collisions may 
become inelastic on account of the ion interaction with the solvent. A natural way to 
include such inelastic effects within the BGK operator is to assume a mismatch between 
the equilbrium current J eq = Y,?^ 1 fi 9 v% and the actual current J carried by the ions. The 
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simplest form for such a mismatch is 

J eq = XJ (47) 

where A is a parameter measuring the degree of inelasticity of the 'effective' ion-ion collisions 
(0 < A < 1 for passive solvents, i.e. momentum absorbing media, while A > 1 denotes active 
solvents, i.e. momentum-imparting media, and A = 1 for standard elastic collisions). 

The relation (|4*7jl is readily implemented by replacing u with \u in the second term within 
parentheses on the right hand side of the eq. (|I0|) . With this modification, the BGK collision 
operator contributes a term — uj(J — J eq ) = —uj{\ — A) J to the momentum equation, so that 
inelastic collisions make themselves felt through an effective frequency u)\ = (1 — X)u. As 
a result, setting aside boundary conditions, i.e. in an infinitely long channel, for simplicity, 
the steady-state current is given by 

= nE/j 
l + w A /7 

This expression has been checked against numerical simulations (with periodic boundary 
conditions) for A = (fully inelastic) and A = 1 (fully elastic), with 7 = 0.1 and uj/j = 
0.1, 1, 10. In Figure 4, we report the ratio J(A = 0)/J(A = 1) as a function of the electric field 
strength, and compare it with the analytical result 7/(7 + uj). From this figure, excellent 
agreement with the analytical results is clearly appreciated. Although the present test 
only serves the purpose of illustrating the basic idea, it is hoped that combination of the 
lattice Fokker-Planck equation with inelastic BGK operators, will permit to explore complex 
situations in one, two and three dimensions, out of reach of analytical methods, such as 
multicomponent fluids (ions, cations and solvent) with heterogeneous surface interactions, 
trapping effects and related phenomena. 

VII. CONCLUSION 

In conclusion, we have developed a lattice version of the Fokker-Planck equation (FPE) 
which, by construction, can solve near-equilibrium kinetic problems in the full single-particle 
phase space. The main scope of the lattice FPE is to account for dissipative interactions 
not resolved at the molecular scale, such as fluid interactions with solid walls and/or solute- 
solvent collisions. The lattice FPE has been applied to the study of electrorheological trans- 
port of a one-dimensional charged fluid, and found to yield satisfactory agreement with a 
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recent non-trivial analytical solution, especially for the five-speed case. In particular, the 
lattice FPE proves capable of predicting the saturation effect resulting from the non-linear 
interaction between the electric field and the constant-flux boundary conditions imposed by 
the presence of equilibrium reservoirs at the channel boundaries. 

The present lattice FPE extends straightforwardly to higher dimensions and it might 
prove useful for the numerical investigation of more complex situations, such as heterogeneus 
channels with kinetic traps, and/or multicomponent fluids, for which analytical solutions are 
no longer available. 
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Appendix A: STABILITY ANALYSIS OF THE LATTICE FPE 

Being derived from the same principles, the lattice Fokker-Planck equation shares many 
features with the Lattice Boltzmann equation. Among the main advantages: i) a very 
efficient sampling of velocity space, which permits to work in steps of size Av = vt rather 
than of a fraction thereof; ii) since the discrete speeds Vi are constant, the streaming operator 
vd x f can be integrated exactly along the characteristics Ax = ViAt, i.e. by a mere shift of 
the discrete distribution from site x to site x + ViAt, iii) the collision operator is completely 
local in space, which makes the Lattice FPE well suited to parallel computing. Of course, 
there are limitations too. In particular, the damping rate is approximately bounded within 
the range 

a/v T < 7 < 1/At 

The lower bound relates to the fact that damping rates below cl/vt imply that particles may 
acquire a drift speed u = a/7 larger than the thermal speed. This endangers the stability 
of the lattice FPE, because the condition u > vt may violate the positive-definiteness of 
the truncated distribution feq IT9|) . which consists only of a very limited number of Hermite 
polynomials. Thus, in full analogy with LBE, the lattice FPE is best suited to describe near- 
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equilibrium, low-Mach number flows, with u/vt -C 1. The upper boundary is dictated by a 
stability condition of the time-marching scheme, specifically of the explicit time integration 
of the collision operator (we recall that the streaming operator is integrated exactly). The 
stability of the lattice FPE can be estimated by means of the following inequality 

,C k At, 

I-tH <! (Al) 
fk 

which states that the change in the macroscopic moment F k (x,t) due to collisions in a time 
step At, should not exceed the value of Fk(x,t) itself. 

In order to elaborate this condition further, we recall the expression of the five Hermite 
polynomials relevant to the Q$ lattice, that is: 

h (y i ) = l i = [1,1, 1,1,1] 

h 1 {v i )=v i = [-2,-l,0,+l,+2] 
h 2 (vi) =v] -v 2 = [3,0,-1,0,3] 
h 3 (vi) = vf - 3v 2 Vl = [-2, 2, 0, -2, 2] 
h A (vi) = vf - 4v 2 v 2 + 4 = [1, -2, 1, -2, 1] 
The corresponding macroscopic moments in equations (J30|) - ([3i|) . are given by: 

F Q = n (A2) 

F\ — J (A3) 

F 2 = nu 2 (A4) 

F 3 = Q-3J4 (A5) 

F 4 = R - Anv 2 T u 2 - 3m4 (A6) 

In view of these expressions, the stability condition (|Aljl is readily recast in a more infor- 
mative form as: 

kjAt < 1, k = 1,4 (A7) 

A more specific result can be obtained by analyzing the dispersion relation associated 
with the lattice FPE. Upon Fourier-transforming the lattice FPE eq. lf^ . fi(x,t) = 
Y,k,u fi(k,uj)e I ( kx ~ ult \ we obtain (J denotes the imaginary unit): 

£ [(e -/(^)At _ 1)8 .. _ C .. At]f . = q (Ag) 
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where CV, is the collision matrix associated with the Fokker-Planck operator. This matrix 
can be computed as follows. Consider the definition of the spectral coefficient C\ in eq. (|23|) . 
and express the function C FP [f] as Cf, where C = d Va [R a + Dd Va ] is the Fokker-Planck 
operator acting upon f(x,v,t). By expanding the function / on the Hermite basis as given 
by eq.fJUU, eq.JSBJ) takes the form 

C k = J2 C kiFi (A9) 
i 

where 

C kl = J h k (v)C wivMv^v (A10) 

is the matrix representation of the Fokker-Planck operator in the global Hermite basis. By 
using eq. (j23J) to express Fi in terms of fj, we obtain: 

Ci = J2 C k ih k {vi)wi h i(vj)fj = C ijfj 

kl j j 

which defines the collision matrix CV, as: 

Cij = WiJ2 h k {vi)C k ihi(vj) (All) 

kl 

The above expression is the operational key of the lattice FPE and lends itself to a fairly 
transparent physical interpretation. The matrix element Cy, expressing the effect of popu- 
lation fi on population fj, is a weighted average over all spectral modes, the weights being 
the Hermite eigenfunctions evaluated at v — v i and v = Vj respectively. 

It is now convenient to normalize the Hermite coefficients h k {yi) as follows: 

h ik = h k {vi)/yjH k 

where H k = J2i h k {vi)wih k {vi) are the normalization factors, explicitly H = 1, Hi = 1, H 2 = 
2, H% = 2, H/i = 2. 

The expressions (jAllj) and (|A10j) are nothing but the matrix representations of the 
Fokker-Planck operator in the local (Dirac's deltas) and global (Hermite polynomials) basis 
functions. The transformation between these two representations is performed by the matrix 
hik- It is readily checked that this matrix fulfills the orthonormality condition 

hi k Wih k j = 5ij 

k 

Consequently, the matrices and C k i are related by a similarity transformation, hence 
they are iso-spectral (they share the same eigenvalues). 
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The spectrum of Cm can be obtained by direct inspection of eq.s IpjUJl - ljHlJ) for the case 
E — 0, and taking into account the definitions (|A2J) . The resulting matrix Cm is identified 
as: 



/(DO 

-7 





-37 
4 7 -4 7 J 
which delivers the following five eigenvalues: 






2 7 









V 



Az 



0,4 



Note the zero eigenvalue associated with mass conservation. 

Acc 0r d ing to staodard a rgum eot S of Lattice Bolt M t h eo ry Q, toe stability coaditioa 
associated with the dispersion relation (|A8j) . reads as follows: 



|1 - kjAt\ < 1, 



1,4 



This is satisfied within the range 



< 7At < 1/2 



(A12) 



(A13) 



It is perhaps interesting to observe that the procedure outlined in this Appendix can 
be applied to a very broad class of kinetic equations, including the Klein-Gordon and 
Schrodinger equations of quantum mechanics [24]. 
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Figure 1: Ratio of numerical to analytical density profiles along a channel of length 10000, for 
a* = 1.0 and a* = 3.0. The solid line and the dashed line correspond to the Q3 and the Q5 models 
respectively. 
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Figure 2: Reduced current J* as a function of the reduced applied field a* = E/Et- Filled symbols 
correspond to the Q3 model while open symbols to the Q5 model. Circles, squares and diamonds 
correspond to 7* = 0.1, 5 and 10 respectively. Solid, dashed and long-dashed lines are the theoretical 
predictions for the three values of 7* . The horizontal line highlights the limiting plateau J* = 1. 
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Figure 3: The pointwise error of the numerical vs. analytical solution for the density and current 
at the midpoint of the channel for the Q$ lattice as a function of the number of grid points N. 
Circles and squares refer to density and current respectively. The main parameters are a* = 1 and 
7* = 10. The dashed lines correpond to N -1 and N~ 1 / 2 convergence. 
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Figure 4: Ratio of currents for elastic vs. inelastic collisions, J\=o/ J\=i, for E 1 /Et = 100. Three 
values of X/lo = 0.1, 1, and 10 correspond to the upper, middle and lower curves respectively. The 
horizontal lines represent the theoretical prediction. 
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